include("PreProcess.jl")
include("EulerRhs.jl")

function main()
    # Order of interpolate polynomials used for approximation
    N = 4

    # Set simulation type
    simtype = "InsentroicVortex"
    if simtype == "InsentroicVortex"
        meshfile = "vortexA04.neu"
    else
        print("Simulation case using unknown")
        exit()
    end

    # Read mesh from meshfile
    Nv, K, Vx, Vy, EToV, BCType = MeshReaderGambit2D(meshfile)

    # Definition of constants
    Nfaces=3; Nfp = N+1; Np = convert(Int, (N+1)*(N+2)/2); NODETOL = 1e-12

    # Compute nodal Set
    x_s,y_s = Nodes(N); r,s = xytors(x_s,y_s)
    V = Vandermonde2D(N, r, s); invV = inv(V)
    MassMatrix = invV'*invV
    Dr,Ds = Dmatrices2D(N, r, s, V)

    # Build coordinates of all the nodes
    va = EToV[:,1]'; vb = EToV[:,2]'; vc = EToV[:,3]'
    x = 0.5*(-(r+s)*Vx[va] + (1 .+ r)*Vx[vb] + (1 .+ s)*Vx[vc])
    y = 0.5*(-(r+s)*Vy[va] + (1 .+ r)*Vy[vb] + (1 .+ s)*Vy[vc])

    # Find all the nodes that lie on each edge
    fmask1 = findall(abs.(s[:] .+ 1) .< NODETOL)
    fmask2 = findall(abs.(r[:] .+ s[:]) .< NODETOL)
    fmask3 = findall(abs.(r[:] .+ 1) .< NODETOL)
    Fmask = [fmask1 fmask2 fmask3]
    Fx = x[Fmask[:],:]; Fy=y[Fmask[:],:]

    # Creat surface integral terms
    LIFT = Lift2D(N,Np,Nfaces,Nfp,Fmask,r,s,V)

    # Calculate geometric factors
    rx,sx,ry,sy,J = GeometricFactors2D(x,y,Dr,Ds)
    nx, ny, sJ = Normals2D(x,y,Dr,Ds,Fmask,Nfp,K)
    Fscale = sJ./J[Fmask[:],:]

    # Build connectivity matrix
    EToE, EToF = tiConnect2D(EToV,Nfaces,K,Nv)
    mapM, mapP, vmapM, vmapP, vmapB, mapB = BuildMaps2D(N,K,Fmask,EToE,EToF,EToV,Vx,Vy,x,y)

    # Compute weak operators (could be done in preprocessing to save time)
    Vr, Vs = GradVandermonde2D(N, r, s)
    Drw = (V*Vr')/(V*V'); Dsw = (V*Vs')/(V*V')
    mapI, mapO, mapW, mapF, mapC, mapD, mapN, mapS = BuildBCMaps2D(Nfp,BCType,vmapM)

    # Initilize the solutions
    Q = InitCondition(Np,K,x,y,0,simtype)

    # Compute initial time-step
    gamma = 1.4
    dt = EulerDt(N,Q,vmapM,gamma,Fscale)
    time = 0; tstep = 1; FinalTime = 1.0

    # Define RK coefficents
    rk4a = [0.0,
            -567301805773.0/1357537059087.0,
            -2404267990393.0/2016746695238.0,
            -3550918686646.0/2091501179385.0,
            -1275806237668.0/842570457699.0];
    rk4b = [1432997174477.0/9575080441755.0,
            5161836677717.0/13612068292357.0,
            1720146321549.0/2090206949498.0,
            3134564353537.0/4481467310338.0,
            2277821191437.0/14882151754819.0]
    rk4c = [0.0,
            1432997174477.0/9575080441755.0,
            2526269341429.0/6820363962896.0,
            2006345519317.0/3224310063776.0,
            2802321613138.0/2924317926251.0];

    # Storage for low storage RK time stepping
    rhsQ = zeros(Float64, Np,K,4)
    resQ = zeros(Float64, Np,K,4)

    # Time loop
    while time < FinalTime
        # Check to see if we need to adjust for final time step
        if time+dt > FinalTime
            dt = FinalTime-time
        end

        # SSK loop
        for Index = 1:5
            # Compute right hand side of compressible Euler equations
            rhsQ  = EulerRhs(Q,N,K,vmapM,vmapP,mapI, mapO, mapW, mapC,Fx,Fy,Drw,Dsw,rx,sx,ry,sy,nx,ny,time,gamma,LIFT,Fscale)

            # Initiate and increment Runge-Kutta residuals
            resQ = rk4a[Index]*resQ + dt*rhsQ

            # Update fields
            Q = Q+rk4b[Index]*resQ
        end

        # Update time and compute new timestep
        time = time+dt; tstep = tstep+1
        dt = EulerDt(N,Q,vmapM,gamma,Fscale)
    end
    Q
end
